SUBROUTINE moments(n, p, fval)
  USE commonvars

  IMPLICIT NONE
  INCLUDE 'mpif.h'

  INTEGER, INTENT(IN)     :: n
  REAL(8), INTENT(IN)     :: p(n)
  REAL(8), INTENT(OUT)    :: fval
  INTEGER                 :: i, tt, ii, Npoly, nprocs_t, myrank_t
  INTEGER, PARAMETER      :: Nfinest = Nfines-1
  REAL(8)                 :: temp(Nmom,Nmomdata), diffmom(Nmom,1), qcons_w_riv, t1, t2, sumbetadf, uu, pfine_temp, pfine_temp1, &
                             cdf_fine(Nfinest), true_cdf_fine(Nfinest), frac_fine_temp, frac_finet(Nfinest), &
                             prob_a_temp, prob_a_temp1, cdf_a(Nab), true_cdf_a(Nab), a_discr_temp, beta, res_prod_func

  EXTERNAL                :: utilfun, emaxfun, simulatemun
  
  CALL cpu_time(t1)

  IF (flag_policy_evaluation == 0 .OR. (flag_policy_evaluation == 1 .AND. flag_optimal_auditprob == 1)) THEN
     nprocs_t = nprocs
     myrank_t = myrank
  ELSE
     nprocs_t = nprocs_new
     myrank_t = my_new_comm_rank
  END IF
  
  iterations=iterations+1
  IF(flag_bootstrap == 0 .AND. myrank_t == 0) write(*,*) 'iteration', iterations

  ii = 0

  ii = ii + 1
  ! theta is the relative value of public consumption
  IF (ii <= nparam) THEN
     IF (flag_sa == 0) THEN
        theta = factor_param(ii)*(ubp(ii)*EXP(p(ii))+lbp(ii))/(1d0+EXP(p(ii)))
     ELSE
        theta = factor_param(ii)*p(ii)
     END IF
  ELSE
     theta = 0.5d0
  END IF

  ii = ii + 1
  !theta_v is the relative value of public consumption for voters; we switched to a model with only one type in terms of relative value
  IF (ii <= nparam) THEN
     IF (flag_sa == 0) THEN
        delta_fine = factor_param(ii)*(ubp(ii)*EXP(p(ii))+lbp(ii))/(1d0+EXP(p(ii)))
     ELSE
        delta_fine = factor_param(ii)*p(ii)
     END IF
  ELSE
     delta_fine = 0.1d0
  END IF

  ii = ii + 1
  IF (ii <= nparam) THEN
     IF (flag_sa == 0) THEN
        cost_run_frac = factor_param(ii)*(ubp(ii)*EXP(p(ii))+lbp(ii))/(1d0+EXP(p(ii)))
     ELSE
        cost_run_frac = factor_param(ii)*p(ii)
     END IF
  ELSE
     cost_run_frac = 0.1d0
  END IF    

  ii = ii + 1
  ! upower
  IF (ii <= nparam) THEN
     IF (flag_sa == 0) THEN
        upower = factor_param(ii)*(ubp(ii)*EXP(p(ii))+lbp(ii))/(1d0+EXP(p(ii)))
     ELSE
        upower = factor_param(ii)*p(ii)
     END IF
  ELSE
     upower = 1d0
  END IF

  ii = ii + 1
  ! Electoral appeal parameter
  IF (ii <= nparam) THEN
     IF (flag_sa == 0) THEN
        eta = factor_param(ii)*(ubp(ii)*EXP(p(ii))+lbp(ii))/(1d0+EXP(p(ii)))
     ELSE
        eta = factor_param(ii)*p(ii)
     END IF
  ELSE
     eta = 1d0
  END IF

  ii = ii + 1
  IF (ii <= nparam) THEN
     IF (flag_sa == 0) THEN
        sigma_run = factor_param(ii)*(ubp(ii)*EXP(p(ii))+lbp(ii))/(1d0+EXP(p(ii)))
     ELSE
        sigma_run = factor_param(ii)*p(ii)
     END IF
  ELSE
     sigma_run = 1d0
  END IF

  ii = ii + 1
  IF (ii <= nparam) THEN
     IF (flag_sa == 0) THEN
        sigma_steal = factor_param(ii)*(ubp(ii)*EXP(p(ii))+lbp(ii))/(1d0+EXP(p(ii)))
     ELSE
        sigma_steal = factor_param(ii)*p(ii)
     END IF
  ELSE
     sigma_steal = 0d0 ! 0.05142432d0
  END IF

  ii = ii + 1
  ! Constant in the learner effect in the production function
  IF (ii <= nparam) THEN
     IF (flag_sa == 0) THEN
        alpha_ie = factor_param(ii)*(ubp(ii)*EXP(p(ii))+lbp(ii))/(1d0+EXP(p(ii)))
     ELSE
        alpha_ie = factor_param(ii)*p(ii)
     END IF
  ELSE
     alpha_ie = 0d0
  END IF

  ii = ii + 1
  ! Effect of being audited and not having stolen on election probabilities parameter
  IF (ii <= nparam) THEN
     IF (flag_sa == 0) THEN
        alpha_nst = factor_param(ii)*(ubp(ii)*EXP(p(ii))+lbp(ii))/(1d0+EXP(p(ii)))
     ELSE
        alpha_nst = factor_param(ii)*p(ii)
     END IF
  ELSE
     alpha_nst = 0d0
  END IF

  ii = ii + 1
  ! Effect of being audited and having stolen a positive amount on election probabilities parameter
  IF (ii <= nparam) THEN
     IF (flag_sa == 0) THEN
        alpha_ea = factor_param(ii)*(ubp(ii)*EXP(p(ii))+lbp(ii))/(1d0+EXP(p(ii)))
     ELSE
        alpha_ea = factor_param(ii)*p(ii)
     END IF
  ELSE
     alpha_ea = 0d0
  END IF

  ii = ii + 1
  ! Effect of being audited and having stolen on election probabilities parameter
  IF (ii <= nparam) THEN
     IF (flag_sa == 0) THEN
        alpha_q = factor_param(ii)*(ubp(ii)*EXP(p(ii))+lbp(ii))/(1d0+EXP(p(ii)))
     ELSE
        alpha_q = factor_param(ii)*p(ii)
     END IF
  ELSE
     alpha_q = 0d0
  END IF
    
  ii = ii + 1
  ! Effect of being audited and having stolen on election probabilities parameter
  IF (ii <= nparam) THEN
     IF (flag_sa == 0) THEN
        alpha_st = factor_param(ii)*(ubp(ii)*EXP(p(ii))+lbp(ii))/(1d0+EXP(p(ii)))
     ELSE
        alpha_st = factor_param(ii)*p(ii)
     END IF
  ELSE
     alpha_st = 0d0
  END IF

  pappeal(1) = 1d0
      
  ii = ii + 1
  ! Effect of being audited and having stolen on election probabilities parameter
  IF (ii <= nparam) THEN
     IF (flag_sa == 0) THEN
        zeta = factor_param(ii)*(ubp(ii)*EXP(p(ii))+lbp(ii))/(1d0+EXP(p(ii)))
     ELSE
        zeta = factor_param(ii)*p(ii)
     END IF
  ELSE
     zeta = 0d0
  END IF
  
  !We use a log-normal distribution for the fine
  !mean of the fine distribution
  mu_fine = 0.0942693d0 + delta_fine
  IF (flag_hfine_policy == 1) THEN
     IF (per_ch_fine > 0d0) THEN
        mu_fine = mu_fine+ABS(mu_fine*per_ch_fine)
     ELSE
        mu_fine = mu_fine-ABS(mu_fine*per_ch_fine)
     END IF
  END IF
  
  !standard deviation of the fine distribution
  sigma_fine = 0.2837051d0

  !The following calculations produce the discrete approximation for the probability of drawing one of the discrite values, frac_fine(:), and one of the corresponding discrete values, pfine(:)
  !Given the parameters of the fine distribution we derive the support points that gives the best approximation (part B in Kennan)
  ! WE USE Nfinest = Nfines-1 to compute the distribution of fines; we then add one more point that corresponds to the case of no conviction with corresponding probability
  DO i = 1, Nfinest
     ! cdf values at the fixed support points
     true_cdf_fine(i) = (2d0*DBLE(i)-1)/(2d0*DBLE(Nfinest))
     !We compute the inverse of the uniform at true_cdf_fine(i); we consider a uniform defined on [0,mu_fine+sigma_fine*E[stealing in data]]
     CALL normal_01_cdf_inv(true_cdf_fine(i),frac_fine_temp)
     frac_finet(i) = EXP(sigma_fine*frac_fine_temp + mu_fine)
  END DO

  !We compute the probabilities at the chosen points that best approximate the log-normal distribution (part A in Kennan)
  CALL normal_01_cdf((LOG(frac_finet(1)) - mu_fine)/sigma_fine,pfine_temp)
  DO i = 2, Nfinest
     CALL normal_01_cdf((LOG(frac_finet(i)) - mu_fine)/sigma_fine,pfine_temp1)
     cdf_fine(i-1) = (pfine_temp + pfine_temp1)/2d0
     pfine_temp = pfine_temp1
  END DO
  cdf_fine(Nfinest) = 1d0

  ! This corresponds to the case of no conviction
  frac_fine(1) = -1d0
  frac_fine(2:Nfines) = frac_finet(:)

  !We compute the probability of the discrete fine values in our grid
  ! pfine(1) corresponds to the probability of no conviction
  pfine(1) = 1d0 - prob_conv
  pfine(2) = cdf_fine(1)
  DO i = 2, Nfinest
     pfine(i+1) = cdf_fine(i) - cdf_fine(i-1)
  END DO
  pfine(2:Nfines) = pfine(2:Nfines)*prob_conv
  
  IF (flag_bootstrap == 0) THEN
     IF(myrank_t == 0) write(*,'(a,3f20.10)') 'cdf_fine', true_cdf_fine(:)
     IF(myrank_t == 0) write(*,'(a)') ''
     IF(myrank_t == 0) write(*,'(a,3f20.10)') 'cdf_fine', cdf_fine(:)
     IF(myrank_t == 0) write(*,'(a,4f20.10)') 'frac_fine', frac_fine(:)
     IF(myrank_t == 0) write(*,'(a,5f20.10)') 'pfine', pfine(:), SUM(pfine(:))    
  END IF
  
  !The ability of mayors is drawn from LOG-N(mu_a,sigma_a)
  !Mean of the ability distribution, we normalize the constant in the production function to 0 and assign the entire constant in the IV estimation to mu_a
  res_prod_func = 3.0858715d0  !3.078571d0
  mu_a = -3.516575d0 + res_prod_func ! -0.2941902d0 !-0.32856629d0 !(ubp(1)*EXP(p(1))+lbp(1))/(1d0+EXP(p(1))) !-.53632455d0
  IF (flag_change_ability_policy == 1) THEN
     IF (ch_mu_a > 0d0) THEN
        mu_a = mu_a + ABS(mu_a*ch_mu_a)
     ELSE
        mu_a = mu_a - ABS(mu_a*ch_mu_a)
     END IF
  END IF
  
  !standard deviation of the ability distribution
  sigma_a = 0.4391229d0  !0.39916454d0 !0.44185212d0 ! .35204222d0 !.3438509d0 !.4730158d0 ! .48566405d0

  IF (flag_change_ability_01stdev == 1) THEN
     mu_a = mu_a + 0.1d0*sigma_a
  END IF
  IF (flag_change_ability_2ndterm_mayors == 1) THEN
     mu_a = -0.401549866285479d0
  END IF

  !Mean of the ability distribution in the second term = mu_a + beta' - beta, where beta' and beta are the estimated constants in the production functions for the first and second term
  beta = mu_a !-.2941902d0 !-0.32856629d0
     
  !parameters of the production function for public consumption                                         
  !we use the same constant in both terms
  !power parameter for public inputs
  alphaq(1) = 0.1941907d0  !.1526144d0 !0.1518577d0
  !power parameter for private inputs
  alphaq(2) = 0.4339804d0  !.429166d0 !0.3725852d0
  !linear term in education
  alphaq(3) = 0d0  !0.0430096d0 !0.0340758d0
  !We normalize the constant to be equal to 0
  alphaq(4) = beta - mu_a
  !flexible polynomial in population
  alphaq(5) = 0d0 !6.66d-06 
  alphaq(6) = 0d0 !-2.26d-10
  alphaq(7) = 0d0 !1.46d-15
  alphaq(8) = 0d0 !-3.75d-21
  alphaq(9) = 0d0 !3.16d-27

  IF (flag_bootstrap == 1) THEN

     alphaq(1) = prod_func_param(1)
     alphaq(2) = prod_func_param(2)
     alphaq(3) = prod_func_param(3)
     alphaq(4) = prod_func_param(4)
     alphaq(5) = prod_func_param(5)
     alphaq(6) = prod_func_param(6)
     alphaq(7) = prod_func_param(7)
     alphaq(8) = prod_func_param(8)
     alphaq(9) = prod_func_param(9)

     mu_a = prod_func_param(10) + res_prod_func
     sigma_a = prod_func_param(11)

!!$     IF (myrank_t==0) write(*,'(a,11f20.10)') 'prod funct param, bootstrap', alphaq(:), mu_a, sigma_a
     
  END IF

  !Polynomial of order n in population
  Npoly = Nparq - 4
  polypop = 0d0
!!$  DO imun = 1, Nmun
!!$     DO i = 1, Npoly
!!$        polypop(imun) = polypop(imun) + alphaq(4+i)*popsize(imun)**DBLE(i)
!!$        !if (myrank_t ==0) write(*,'(a,2i4,4f20.10)') 'polypop', imun, i, polypop(imun), alphaq(4+i), popsize(imun), popsize(imun)**DBLE(i)
!!$     END DO
!!$  END DO
  
  !Ea and Eap are the expectation of ability in the first and second term (ability is log-normally distributed
!!$  Ea = EXP(mu_a+sigma_a**2d0/2d0)
!!$  Eap = EXP(mu_ap+sigma_ap**2d0/2d0)

!!$. ivregress gmm lngdp_pc_01_04 lindustry_index pref_edu PesResAlfaRate* pib2001* d_2nd_term (lpinputs_pc_01_04 = fpm_d1_200* fpm_d2_200* fpm_d3_200*  fpm_d4_200*  fpm_d5_200* fpm_d6_200* fpm_d7_200* fpm_d8_200* fpm_d9_200* fpm_d10_200* fpm_d11_200* fpm_d12_200* fpm_d1_pop_200* employment_rate2000 wage_priv nfirms ln_pop_01_04) if nterms >= 1, rob
!!$
!!$Instrumental variables (GMM) regression           Number of obs   =        475
!!$                                                  Wald chi2(8)    =    8471.29
!!$                                                  Prob > chi2     =     0.0000
!!$                                                  R-squared       =     0.6210
!!$GMM weight matrix: Robust                         Root MSE        =     .45851
!!$
!!$------------------------------------------------------------------------------------
!!$                   |               Robust
!!$    lngdp_pc_01_04 |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
!!$-------------------+----------------------------------------------------------------
!!$ lpinputs_pc_01_04 |   .1343342   .0341213     3.94   0.000     .0674577    .2012107
!!$   lindustry_index |   .4431576   .0923581     4.80   0.000     .2621389    .6241762
!!$          pref_edu |   .0517908   .0365824     1.42   0.157    -.0199093    .1234909
!!$    PesResAlfaRate |   .0003895   .0000202    19.26   0.000     .0003498    .0004291
!!$PesResAlfaRate_2nd |  -.0000521   .0000228    -2.28   0.022    -.0000968   -7.39e-06
!!$           pib2001 |   2.32e-07   9.47e-08     2.45   0.014     4.67e-08    4.18e-07
!!$       pib2001_2nd |   4.83e-07   1.75e-07     2.76   0.006     1.40e-07    8.27e-07
!!$        d_2nd_term |   .3850212   .1656483     2.32   0.020     .0603565    .7096859
!!$             _cons |  -3.441379   .2032496   -16.93   0.000    -3.839741   -3.043017
!!$------------------------------------------------------------------------------------

!!$. scalar beta_1 = const + Pes*E_Pes + pib*E_pib
!!$
!!$. scalar beta_2 = const + d_2nd + Pes*E_Pes + Pes_2nd * E_Pes_2nd + pib*E_pib + pib_2nd*E_pib_2nd
!!$
!!$. di beta_1
!!$-.35443026
!!$
!!$. di beta_2
!!$-.32818725

    !The following calculations produce the discrete approximation for the probability of drawing one of the discrite values, prob_a(:), and one of the corresponding discrete values, a_discr(:)
  !Given the parameters of the ability distribution we derive the support points that gives the best approximation (part B in Kennan)
  DO i = 1, Nab
     ! cdf values at the fixed support points
     true_cdf_a(i) = (2d0*DBLE(i)-1)/(2d0*DBLE(Nab))
     !We compute the inverse of the standard normal at true_cdf_a(i)
     CALL normal_01_cdf_inv(true_cdf_a(i), a_discr_temp)
     a_discr(i) = sigma_a*a_discr_temp + mu_a
  END DO
  
  !We compute the probabilities at the chosen points that best approximate the log-normal distribution (part A in Kennan)
  CALL normal_01_cdf((a_discr(1) - mu_a)/sigma_a,prob_a_temp)
  DO i = 2, Nab
     CALL normal_01_cdf((a_discr(i) - mu_a)/sigma_a,prob_a_temp1)
     cdf_a(i-1) = (prob_a_temp + prob_a_temp1)/2d0
     prob_a_temp = prob_a_temp1
  END DO
  cdf_a(Nab) = 1d0

  !We compute the probability of the discrete fine values in our grid
  pability(1) = cdf_a(1)
  DO i = 2, Nab
     pability(i) = cdf_a(i) - cdf_a(i-1)
  END DO

  !We take the exp of a_discr(i) to make it a log-normal
  a_discr(:) = EXP(a_discr(:))

  IF (flag_bootstrap == 0) THEN
     IF(myrank_t == 0) write(*,'(a)') ''
     IF(myrank_t == 0) write(*,'(a,3f20.10)') 'true cdf ability', true_cdf_a(:)
     IF(myrank_t == 0) write(*,'(a,3f20.10)') 'discrete ability', a_discr(:)
     IF(myrank_t == 0) write(*,'(a,3f20.10)') 'cdf ability', cdf_a(:)
     IF(myrank_t == 0) write(*,'(a,3f20.10)') 'prob discrete ability', pability(:)
  END IF
  
  !parameters of the wage function for past mayors
  !constant
  alphaw(1) = -3.086317d0
  !education
  alphaw(2) = 4.096486d0
  !age
  alphaw(3) = 1.035222d0
  !age^2
  alphaw(4) = -.0975276d0
  !audit*ln(popsize)
  alphaw(5) = 0d0
  !audit*d_steal*ln(popsize)
  alphaw(6) = 0d0
  !dummy for medium size munic
  alphaw(7) = 1.144047d0
  !dummy for large munic
  alphaw(8) = 3.215137d0
  !standard deviation
  sigmaw = 6.335244d0 !6.464272d0

!!$ reg wage_mean_past_norm  pref_edu age age2 pop_d2 pop_d3 if wage_mean_past_norm >.5668408  & wage_mean_past_norm <  10000 & age<99, rob
!!$
!!$Linear regression                               Number of obs     =      3,725
!!$                                                F(5, 3719)        =     143.06
!!$                                                Prob > F          =     0.0000
!!$                                                R-squared         =     0.1392
!!$                                                Root MSE          =     6.3404
!!$
!!$------------------------------------------------------------------------------
!!$             |               Robust
!!$wage_mean_~m |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
!!$-------------+----------------------------------------------------------------
!!$    pref_edu |   4.096486   .1794472    22.83   0.000     3.744662    4.448311
!!$         age |   1.035222   .1830252     5.66   0.000     .6763825    1.394062
!!$        age2 |  -.0975276   .0182942    -5.33   0.000    -.1333952     -.06166
!!$      pop_d2 |   1.144047    .211215     5.42   0.000     .7299387    1.558156
!!$      pop_d3 |   3.215137   .3393853     9.47   0.000     2.549737    3.880536
!!$       _cons |  -3.086317   .5092068    -6.06   0.000    -4.084668   -2.087965
!!$------------------------------------------------------------------------------

  !Standard deviation of the private input shocks
  sigma_zzpr = 0d0 !(ubp(20)*EXP(p(20))+lbp(20))/(1d0+EXP(p(20)))

!!$  !Standard deviation of the shocks to the decision to run
!!$  sigma_run = 1d0 ! 1d0 !0.1d0 ! (ubp(19)*EXP(p(19))+lbp(19))/(1d0+EXP(p(19)))
  
  ! WE WILL ASSUME LOG PREFERNCES
  !preference parameters                                                                                
  delta = 0d0 ! 2d0 ! 0d0 ! p(6) if we use log utility for public consumption, otherwise we have to estimate the parameter 
  !gamma = 1 if we use log utility for private consumption, otherwise we have to estimate the parameter 
  gamma = 2d0

  ! WE SHOULD GET RID OF THE MEASURMENT ERRORS
  !Mean of the measurement error in stealing; a positive mean means that it is more likely that a positive clearical mistake (leading to corruption) is made than negative (leading to negative corruption)
  !As a percentage of funds (see simulatemun)
  mu_steal = 0d0 !(ubp(22)*EXP(p(22))+lbp(22))/(1d0+EXP(p(22)))

  !We transform the level of public consumption to take into account the degree of rivalry
  !This is also the variable we output for public consumption
  !If we divide by (popsizeg)**eta, the utility at subsistence level goes to -infinity 
  qcons_w_riv = subqcons_vec(2) + small_no

  subutil = 0d0
  CALL utilfun(subpcons,qcons_w_riv,theta,uu)
  IF (flag_bootstrap == 0 .AND. myrank_t==0) WRITE(*,*) 'uu,subpcons,subqcons,theta', uu,subpcons,qcons_w_riv,theta
  sumbetadf = 0d0
  DO tt = 1, Tend
     sumbetadf = sumbetadf + betadf**(tt-1)
     subutil(Tend - tt + 1) = uu*sumbetadf
  END DO

!!$  subutil = (/ -9146.35964062422, -8312.61187986127, -7461.84885867459, -6593.72332685144, -5707.88094744007, -4803.96015212234, -3881.59199363487, -2940.39999517826, -1979.99999675314, -999.999998360174 /)

  IF (flag_bootstrap == 0 .AND. myrank_t==0) WRITE(*,*) 'subutil', subutil(:)

  IF (flag_bootstrap == 0 .AND. myrank_t == 0) THEN
     write(*,*) ''
     write(*,*) 'param', p(:)
     write(*,*) 'electoral probability: alpha_ea, alpha_st, alph_q, alpha_ie, alpha_nst', alpha_ea, alpha_st, alpha_q, alpha_ie, alpha_nst
     write(*,*) 'pappeal', pappeal
     write(*,*) 'mu_fine, sigma_fine, delta_fine', mu_fine, sigma_fine, delta_fine
     write(*,*) 'plability', pability
     write(*,*) 'prod fn: pu*ab, pr*ab, educ, const', alphaq(1:4)
     write(*,*) 'prod fn: poly pop', alphaq(5:Nparq)
     write(*,*) 'theta', theta
     write(*,*) 'upower', upower
     write(*,*) 'wages past mayors', alphaw
     write(*,*) 'delta, gamma', delta, gamma
     write(*,*) 'sigma_run, sigma_zzpr', sigma_run, sigma_zzpr
     write(*,*) 'sigma_steal, mu_steal', sigma_steal, mu_steal
     write(*,*) 'cost_run_frac, upower', cost_run_frac, upower
     write(*,*) 'mu_a, sigma_a', mu_a, sigma_a
     write(*,*) 'beta', beta
     write(*,*) 'eta', eta
     write(*,'(a,11f18.14)') 'param inputs', theta, delta_fine, cost_run_frac, upower, eta, sigma_run, sigma_steal, alpha_ie, alpha_nst, alpha_ea, alpha_q
  END IF

  CALL emaxfun()
  
  CALL simulatemun()

  !----------- COMPUTE CRITERION FUNCTION ------------!

  ! Compute the difference in moments
  diffmom(1:Nmom,1) = simmom(1:Nmom,1) - datamom(1:Nmom,1)

  DO i=1,Nmom
     temp(i,1)=SUM(diffmom(:,1)*INV_V(:,i))
     if (flag_bootstrap == 0 .AND. myrank_t == 0) write(*,'(A,1I3,x,A,1F12.6,x,A,1F12.6,x,A,1F12.6,x,A,1F12.6,x,A,1F12.6,x,A,1F12.6)') 'i',i,'dfval',diffmom(i,1)*INV_V(i,i)*diffmom(i,1),'simmom',simmom(i,1),'datamom',datamom(i,1),'diff',diffmom(i,1),'V',V(i,i),'weight',INV_V(i,i)
  END DO
  if (flag_bootstrap == 0 .AND. myrank_t == 0) write(*,*) ''

  fval=SUM(temp(:,1)*diffmom(:,1))

  CALL cpu_time(t2)
  IF (flag_bootstrap == 0 .AND. myrank_t == 0) THEN
     write(*,*) 'simmom', simmom
     write(*,*) 'datamom', datamom
     write(*,*) 'fval', fval
     write(*,*) 'time for one set of parameters', t2-t1
  END IF

END SUBROUTINE moments
